Energy Burden

What is energy burden?

Energy burden refers to the percent of a household’s income spent on energy costs.

An individual household’s energy burden might be high if they have low income, meaning the electric bill eats up more of each paycheck than a high income household. Or, someone’s energy burden might higher if they live poorly insulated house, their heating/cooling system is inefficient or outdated, they live in a climate that demands high heating or cooling, or all of the above.

The median energy burden for The United States is 3.1%.¹ The average energy burden for the state of Virginia is slightly lower at 2%.²

Average energy burden for each census tract in the Charlottesville region ranges from about 1.5% to 6.5%. However, averaging all households in a tract makes it harder to see variation within one tract. So instead, we use 6% as a cutoff to define energy burdened households, and then we calculate that percentage of burdened households in each tract. The map below shows the percentage of households in each census tract that are energy burdened (meaning they spend 6% or more of their income on energy costs).

pal <- colorNumeric("Blues", reverse = FALSE, domain = cvlshapes$percentburdened)
leaflet(cvlshapes) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(data = cvlshapes,
              fillColor = ~pal(percentburdened),
              weight = 1,
              opacity = 1,
              color = "white",
              fillOpacity = 0.6,
              highlight = highlightOptions(weight = 2, fillOpacity = 0.8, bringToFront = T),
              popup = paste0(cvlshapes$NAME.y, "<br>",
                            "Pct. Burdened: ", round(cvlshapes$percentburdened, 2))) %>%
  addLegend("bottomright", pal = pal, values = cvlshapes$percentburdened,
            title = "Percent of Energy<br>Burdened Households", opacity = 0.7)

As we can see from this map, the places with the highest number of energy burdened households are in the city of Charlottesville and the western parts of Albemarle County and Louisa County. There is also a huge amount of variation across the region.

What factors are related to energy burden?

Now that we know what energy burden means and what it looks like in the Charlottesville region, the next question is what demographic factors influence energy burden. To find this out, we took the variables below and created terciles, where the observations are arranged from lowest to highest and divided equally into three, ranked groups. 1 indicates the lowest values, and 3 indicates the highest. We then calculated the average percentage of energy burdened households in each group.

Socio-economic factors

The following graphs shows the association between poverty, unemployment, and median household income with energy burden. In each graph, a low rank (1) indicates the tracts with the lowest rate of each variable and high (3) indicates those with the highest.

As we can see from these graphs, higher poverty rates, higher unemployment rates, and lower median household income are all associated with a higher percentage of energy-burdened households.

dat <- all

## poverty plot 
dat1 <- bi_class(dat, x = povrateE, y = percentburdened, style = "quantile", dim = 3)
dat1$povrank <- stri_extract(dat1$bi_class, regex = '^\\d{1}(?=-\\d)')

dat2 <- dat1 %>%  
  dplyr::select(allhseE, percentburdened, povrank) %>%
  filter(!is.na(percentburdened), !is.na(povrank)) %>% 
  mutate(mh_num = allhseE*(percentburdened/100)) %>% 
  group_by(povrank) %>% 
  summarize(mh_num = sum(mh_num),
            num = sum(allhseE),
            mh_per = round((mh_num/num)*100,1))

povplot <- dat2 %>% 
  ggplot(aes(povrank, mh_per, fill=povrank)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(labels = c("low", "mid", "high"),
                    values = c("#a5add3", "#dfd0d6", "#5698b9"),
                    name = "Rank") + 
  labs(x = "% in poverty", y = "% of burdened households") + theme_bw()


## unemployment plot
dat6 <- bi_class(dat, x = unempE, y = percentburdened, style = "quantile", dim = 3)
dat6$unemployrank <- stri_extract(dat6$bi_class, regex = '^\\d{1}(?=-\\d)')

dat6 <- dat6 %>%  
  dplyr::select(allhseE, percentburdened, unemployrank) %>%
  filter(!is.na(percentburdened), !is.na(unemployrank)) %>% 
  mutate(ph_num = allhseE*(percentburdened/100)) %>% 
  group_by(unemployrank) %>% 
  summarize(ph_num = sum(ph_num),
            num = sum(allhseE),
            ph_per = round((ph_num/num)*100,1))

employplot <- dat6 %>% 
  ggplot(aes(unemployrank, ph_per, fill=unemployrank)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(labels = c("low", "mid", "high"),
                    values = c("#a5add3", "#dfd0d6", "#5698b9"),
                    name = "Rank") + 
  labs(x = "% unemployed", y = "% of burdened households") + theme_bw()


## income plot
dat7 <- bi_class(dat, x = hhincE, y = percentburdened, style = "quantile", dim = 3)
dat7$incrank <- stri_extract(dat7$bi_class, regex = '^\\d{1}(?=-\\d)')

dat7 <- dat7 %>%  
  dplyr::select(hhincE, percentburdened, incrank) %>%
  filter(!is.na(percentburdened), !is.na(incrank)) %>% 
  mutate(ph_num = hhincE*(percentburdened/100)) %>% 
  group_by(incrank) %>% 
  summarize(in_num = sum(ph_num),
            num = sum(hhincE),
            in_per = round((ph_num/num)*100,1))

incomeplot <- dat7 %>% 
  ggplot(aes(incrank, in_per, fill=incrank)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(labels = c("low", "mid", "high"),
                    values = c("#a5add3", "#dfd0d6", "#5698b9"),
                    name = "Median household income rank") + 
  labs(x = "Median household income", y = "% of burdened households") + theme_bw()

ggarrange(povplot, employplot, incomeplot, ncol = 3, common.legend = TRUE)

Racial demographics

The following graphs shows the association between race and energy burden. Just like in the graphs above, in each graph, a low rank (1) indicates the tracts with the lowest rate of each variable and high (3) indicates those with the highest.

The graphs show that percent of energy burdened households is much higher in tracts with higher percent non-White (particularly Black) residents, but much lower in tracts with higher percent White populations.

all <- all %>% 
  mutate(bipoc = 100 - whiteE)

# percent non-white plot
dat8 <- bi_class(all, x = bipoc, y = percentburdened, style = "quantile", dim = 3)
dat8$bipocrank <- stri_extract(dat8$bi_class, regex = '^\\d{1}(?=-\\d)')

dat8 <- dat8 %>%  
  dplyr::select(hhincE, percentburdened, bipocrank) %>%
  filter(!is.na(percentburdened), !is.na(bipocrank)) %>% 
  mutate(ph_num = hhincE*(percentburdened/100)) %>% 
  group_by(bipocrank) %>% 
  summarize(in_num = sum(ph_num),
            num = sum(hhincE),
            in_per = round((ph_num/num)*100,1))

bipocplot <- dat8 %>% 
  ggplot(aes(bipocrank, in_per, fill=bipocrank)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(labels = c("low", "mid", "high"),
                    values = c("#a5add3", "#dfd0d6", "#5698b9"),
                    name = "Rank") + 
  labs(x = "% Non-White", y = "% of burdened households") + theme_bw()

# percent Black plot 
dat4 <- bi_class(dat, x = blackE, y = percentburdened, style = "quantile", dim = 3)
dat4$perBlackrank <- stri_extract(dat4$bi_class, regex = '^\\d{1}(?=-\\d)')

dat4 <- dat4 %>%  
  dplyr::select(allhseE, percentburdened, perBlackrank) %>%
  filter(!is.na(percentburdened), !is.na(perBlackrank)) %>% 
  mutate(ca_num = allhseE*(percentburdened/100)) %>% 
  group_by(perBlackrank) %>% 
  summarize(ca_num = sum(ca_num),
            num = sum(allhseE),
            ca_per = round((ca_num/num)*100,1))

raceplot <- dat4 %>%  
  ggplot(aes(perBlackrank, ca_per, fill=perBlackrank)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(labels = c("low", "mid", "high"),
                    values = c("#a5add3", "#dfd0d6", "#5698b9"),
                    name = "Rank") + 
  labs(x = "% Black", y = "% of burdened households") + theme_bw()

# percent White plot
dat9 <- bi_class(dat, x = whiteE, y = percentburdened, style = "quantile", dim = 3)
dat9$perWhiterank <- stri_extract(dat9$bi_class, regex = '^\\d{1}(?=-\\d)')

dat9 <- dat9 %>%  
  dplyr::select(allhseE, percentburdened, perWhiterank) %>%
  filter(!is.na(percentburdened), !is.na(perWhiterank)) %>% 
  mutate(ca_num = allhseE*(percentburdened/100)) %>% 
  group_by(perWhiterank) %>% 
  summarize(ca_num = sum(ca_num),
            num = sum(allhseE),
            ca_per = round((ca_num/num)*100,1))

whiteplot <- dat9 %>%  
  ggplot(aes(perWhiterank, ca_per, fill=perWhiterank)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(labels = c("low", "mid", "high"),
                    values = c("#a5add3", "#dfd0d6", "#5698b9"),
                    name = "Rank") + 
  labs(x = "% White", y = "% of burdened households") + theme_bw()

ggarrange(bipocplot, raceplot, whiteplot, ncol = 3, common.legend = TRUE)

Home ownership

Homeowners have much more incentive to update their to invest in updates to make their homes more energy efficient than are landlords; homeowners can save themselves money on their electric bills whereas landlords would likely only be saving their tenants money by reducing electricity costs. The effect of this difference is obvious when we look at the percentage of energy-burdened renters versus energy-burdened owners. In the graph below, we can see that home ownership actually is correlated with energy burden—in general, renters face far greater energy burdens than owners.

# filter to just variables of interest
bar <- all %>% 
  group_by(county) %>% 
  summarize(pctrenters = mean(percent_burdened_renters),
            pctowners = mean(percent_burdened_owners))
bar <- bar %>% arrange(pctrenters)

# convert from wide to long
bar_long <- gather(bar, type, percent, pctrenters:pctowners)

bar_long$county <- ifelse(bar_long$county == "Charlottesville city", "Charlottesville", bar_long$county)

# graph
ggplot(bar_long, aes(x = county, y = percent, fill = type)) +
  geom_col(position = "dodge") +
  scale_fill_manual(labels = c("% burdened (owners)", "% burdened (renters)"),
                    values = c("#a5add3", "#5698b9"),
                    name = "Resident type") + 
  labs(x = "Locality", y = "% of burdened households") + theme_bw()

Temperature

This map shows the correlation between the number of energy-burdened households and average maximum July temperatures. Tracts are divided into groups based on their relative ranking on percent energy burdened households and the average maximum July temperatures. Mapping both variables together allows the viewer to see how tracts that are hotter, on average, also tend to have a higher percentage of energy-burdened households. This trend is particularly stark in Louisa County.

# manually define palette
bipal <- c("#e8e8e8", "#dfd0d6", "#be64ac", # A-1, A-2, A-3,
           "#ace4e4", "#a5add3", "#8c62aa", # B-1, B-2, B-3
           "#5ac8c8", "#5698b9", "#3b4994") # C-1, C-2, C-2

cvlshapes <- cvlshapes %>%
  mutate(burdenrank = ntile(percentburdened, 3),
         temprank = ntile(July_AvgMaxTF, 3)) %>%
  mutate(burdenrank = if_else(burdenrank == 1, 'A',
                          if_else(burdenrank == 2, 'B', 'C')),
         biclass = paste0(burdenrank, temprank))

# make legend
bipal3 <- tibble(
  "3-3" = "#3b4994", # high average burden, high temperature
  "2-3" = "#8c62aa",
  "1-3" = "#be64ac", # low average burden, high temp
  "3-2" = "#5698b9",
  "2-2" = "#a5add3", # medium avg burden, medium temp
  "1-2" = "#dfd0d6",
  "3-1" = "#5ac8c8", # high avg burden, low temp
  "2-1" = "#ace4e4",
  "1-1" = "#e8e8e8" # low avg burden, low temp
) %>%
  gather("group", "fill")

bipal3 <- bipal3 %>%
  separate(group, into = c("burdenrank", "temprank"), sep = "-") %>%
  mutate(burdenrank = as.integer(burdenrank),
         povertyrank = as.integer(temprank))

legend2 <- ggplot() +
  geom_tile(
    data = bipal3,
    mapping = aes(
      x = burdenrank,
      y = temprank,
      fill = fill)
  ) +
  scale_fill_identity() +
  labs(x = expression("% energy burdened" %->%  ""),
       y = expression("Higher temperature" %->% "")) +
  theme_void() +
  # make font small enough
  theme(
    axis.title = element_text(size = 6),
    axis.title.y = element_text(angle = 90)
  ) +
  # quadratic tiles
  coord_fixed()

# Jacob's fix for the legend
ggsave(plot = legend2, filename = "bivariate_legend.svg",
       width = 1, height = 1)
addLocalLogo <- function(map,
                    img,
                    alpha = 1,
                    src = c("remote", "local"),
                    url,
                    position = c("topleft", "topright",
                                 "bottomleft", "bottomright"),
                    offset.x = 50,
                    offset.y = 13,
                    width = 60,
                    height = 60) {

  if (inherits(map, "mapview")) map <- mapview2leaflet(map)
  stopifnot(inherits(map, c("leaflet", "leaflet_proxy")))

  position <- position[1]
  src <- src[1]
  
  div_topleft <- paste0("newDiv.css({
                        'position': 'absolute',
                        'top': '", offset.y, "px',
                        'left': '", offset.x, "px',
                        'background-color': 'transparent',
                        'border': '0px solid black',
                        'width': '", width, "px',
                        'height': '", height, "px',
});")
  
  div_topright <- paste0("newDiv.css({
                         'position': 'absolute',
                         'top': '", offset.y, "px',
                         'right': '", offset.x, "px',
                         'background-color': 'transparent',
                         'border': '0px solid black',
                         'width': '", width, "px',
                         'height': '", height, "px',
                         });")
  
  div_bottomleft <- paste0("newDiv.css({
                           'position': 'absolute',
                           'bottom': '", offset.y, "px',
                           'left': '", offset.x, "px',
                           'background-color': 'transparent',
                           'border': '0px solid black',
                           'width': '", width, "px',
                           'height': '", height, "px',
                           });")
  
  div_bottomright <- paste0("newDiv.css({
                            'position': 'absolute',
                            'bottom': '", offset.y, "px',
                            'right': '", offset.x, "px',
                            'background-color': 'transparent',
                            'border': '0px solid black',
                            'width': '", width, "px',
                            'height': '", height, "px',
                            });")
  
  div <- switch(position,
                topleft = div_topleft,
                topright = div_topright,
                bottomleft = div_bottomleft,
                bottomright = div_bottomright)
  
  div_funk <- paste0("function(el, x, data) {
                     // we need a new div element because we have to handle
                     // the mouseover output seperately
                     // debugger;
                     function addElement () {
                     // generate new div Element
                     var newDiv = $(document.createElement('div'));
                     // append at end of leaflet htmlwidget container
                     $(el).append(newDiv);
                     //provide ID and style
                     newDiv.addClass('logo');\n",
                     div,
                     "return newDiv;
                     }")
  
  div_add <- paste0("// check for already existing logo class to not duplicate
                    var logo = $(el).find('.logo');
                    if(!logo.length) {
                    logo = addElement();")

  style <- paste0(', style="opacity:',
                  alpha,
                  ';filter:alpha(opacity=',
                  alpha * 100, ');"')
  
  div_html <- paste0("logo.html('<img src=", '"', img, '"',
                     ", width=", width, ", height=", height, style,
                     ", ></a>');
                       var map = HTMLWidgets.find('#' + el.id).getMap();
                     };
                     }")
  
  render_stuff <- paste0(div_funk, div_add, div_html)
  
  map <- htmlwidgets::onRender(map, render_stuff)
  
  map
}

df_4326 <- st_transform(cvlshapes, 4326)

factpal <- colorFactor(bipal, domain = df_4326$biclass)

leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(data = df_4326,
              fillColor = ~factpal(biclass),
              weight = 1,
              opacity = 1,
              color = "white", 
              fillOpacity = 0.8,
              highlight = highlightOptions(
                weight = 2,
                fillOpacity = 0.8,
                bringToFront = T),
              popup = paste0("County: ", df_4326$county, "<br>",
                             "Percent Burdened: ", 
                             round(df_4326$percentburdened, 2), "<br>",
                             "Average Max July Temperature: ", 
                             round(df_4326$July_AvgMaxTF,1))) %>% 
    addLocalLogo("bivariate_legend.svg", src = "local",
          position = "topright", width = 100, height = 100,
          alpha = 0.8)

Click here info on the data sources presented in this story

• Low-Income Energy Affordability Data, 2018 Update
• Daymet: Monthly Climate Summaries on a 1-km Grid for North America
• Census Bureau’s American Community Survey 5-Year Estimates from 2015-2019

---

1. National energy burden average: The American Council for an Energy-Efficiient Economy↩︎

2. Virginia energy burden average: The Office of Energy Efficiency & Renewable Energy↩︎